// Author: Mark Richardson
// Date: 01/26/2022
// Revised date: 08/22/2022 - to make standard deviation on naive prior a parameter passed from R
// Revised date: 10/16/2022 - hierarchical rater model and allow separate traites by party
// Multi-rater item response model for latent traits with naive prior on latent latent traits


data {
  
  int<lower = 0> N; // number of ratings
  
  int<lower = 0> A_fixed; // number of agencies with single trait
  
  int<lower = 0> A_pid; // number of agencies with one trait per party

  int<lower = 0> R; // number of raters
  
  int<lower = 0> D; // number of departments
  
  int<lower = 1, upper = A_pid> aa[N]; // agency index
  
  int<lower = 1, upper = R> rr[N]; // rater index
  
  int<lower = 1, upper = D> dd[R]; // department index
  
  vector[N] y; // vector of ratings
  
  int fixed[N]; // Indicator variable for fixed or partisan
  
  int dem[N]; // Indicator variable for dems

  real<lower = 0> sd_naive; // Naive prior standard deviation
  
}


parameters {
  
  vector[D] mu_alpha; // dept mean for rater intercepts
  
  vector<lower = 0>[D] mu_beta; // dept mean for rater slopes
  
  vector<lower = 0>[D] sigma_alpha; // scale of rater intercepts
  
  vector<lower = 0>[D] sigma_beta; // scale of rater slopes

  vector[A_fixed] x_fixed; // latent performance
  
  vector[A_pid] x_dem; // latent performance - dems
  
  vector[A_pid] x_rep; // latent performance - reps

  real<lower = 0> sigma; // scale of latent performance
  
  //real<lower = 0> sigma_dem; // scale of latent performance

  //real<lower = 0> sigma_rep; // scale of latent performance

  vector[R] alpha_tilde;
  
  vector<lower = 0>[R] beta_tilde;
  
}

transformed parameters {
  
  vector[R] alpha; // rater intercept parameters
  
  vector<lower = 0>[R] beta; // rater slope parameters
  
  for (r in 1:R) {
    
    alpha[r] = mu_alpha[dd[r]] + sigma_alpha[dd[r]] * alpha_tilde[r];
    
   beta[r] = mu_beta[dd[r]] + sigma_beta[dd[r]] * beta_tilde[r];
    
  }
}

model {
  
  mu_alpha ~ normal(0, 5);
  
  mu_beta ~ normal(0, 5);
  
  x_fixed ~ normal(0, sd_naive); // Naive prior on latent trait
  
  x_dem ~ normal(0, sd_naive); // Naive prior on latent trait

  x_rep ~ normal(0, sd_naive); // Naive prior on latent trait
  
  sigma ~ normal(0, 5);
  
  sigma_alpha ~ normal(0, 5);
    
  sigma_beta ~ normal(0, 5);
  
  alpha_tilde ~ std_normal();
  
  beta_tilde ~ std_normal();
  
  for (n in 1:N) {
    
    if (fixed[n] == 1)
      
      y[n] ~ normal(alpha[rr[n]] + beta[rr[n]] * x_fixed[aa[n]], sigma);
      
    else if (dem[n] == 1)
    
      y[n] ~ normal(alpha[rr[n]] + beta[rr[n]] * x_dem[aa[n]], sigma);
      
    else 
      
      y[n] ~ normal(alpha[rr[n]] + beta[rr[n]] * x_rep[aa[n]], sigma);

  }
}

